%{
                                                      analyzeNuclei.m V1
                                                     Written by Rajeev Kant
                                                        Sept. 4th, 2020
                                    Coulombe Lab for Cardiovascular Regenerative Engineering
                            Center for Biomedical Engineering, Brown Univeristy, Providence, RI, 02903

This function takes in a defined ARA region and performs the following
tasks:
1. Perform watershed segmentation on region subROIs
2. Evaluate nuclei based on count and area
3. Calculates nearest neighbors within a specified radius for subROIs
4. Calculates nuclear aspect ratio and filters out non-aligned nuclei
5. Creates summary stats for output and saves appropriate figures/plots

The user may need to edit the following lines: 44 (sharpness), 83 (particle
removal algorithm), 174-179 (nearest neighbor search radius), 212-217
(histology aspect ratio thresholds)

Inputs:
1. region: BIF or BR image
2. ROIpos: ROI position matrix from getROIs.m
3. radRing: Aortic ring radius from ARA_Analyze.m
4. scale: Scale factor
5. areaVal: px value for additional areaopening (if needed)
6. regName: Character array for region ('IF' or 'RM') for naming
7. path: file path
8. filename: name of original input RGB image from ARA_Analyze.m

Outputs:
1. Stats1-6: regionprops stats files for each subROI
2. areaVal: value used in areaopening to apply to next region
3. binNuc: Binary nuclei images for each subROI
4. summary: summary stats for count, nearest 


Original publication for citation: 
Kant, R. J., Bare, C. F., and Coulombe, K. L. K. "Tissues with 
patterned vessels or protein release induce vascular chemotaxis
in an in vitro platform", Tissue Eng. Part A, 2021.
doi:10.1089/ten.TEA.2020.0269

%}
function [stats1, stats2, stats3, stats4, stats5, stats6, areaVal, binNuc, summary] = analyzeNuclei(region, ROIpos, radRing, scale, areaVal, regName, path, filename)
%% Sharpen image to increase contrast

% Prompts are silenced but can be re-enabled if needed
sharp = imsharpen(region,'Radius',3,'Amount',2);
%imshowpair(region, Bsh, 'montage');
%title('Original image on left, sharpened image on right');
%strs  = {'Yes','No'};
%prmpt = 'Use sharpened image?';
%[ans] = listdlg('PromptString',prmpt,'SelectionMode','single','ListString',strs,'ListSize',[120,40]);

%if ans == 1 % Assign sharped image to variable region, otherwise continue with original imported image
    region = sharp;
%end
%close all;

%% Set subROIs

% Set predefined ROI positions based on input image and ROIpos
subR1 = imcrop(region,'Position', ROIpos(1,:));
subR2 = imcrop(region,'Position', ROIpos(2,:));
subR3 = imcrop(region,'Position', ROIpos(3,:));
subR4 = imcrop(region,'Position', ROIpos(4,:));
subR5 = imcrop(region,'Position', ROIpos(5,:));
subR6 = imcrop(region,'Position', ROIpos(6,:));
ROIs = cat(3, subR1, subR2, subR3, subR4, subR5, subR6);
subROI_px = ROIpos(1,3)^2; % subROI area size - subR# size is actually 1 row and column greater than specified ROIpos

%% Binarize and remove small particles based on size exclusion 
% Alter bwareaopen and numel() multiplier as needed

% Initialize variables
segs = zeros(size(ROIs,1),size(ROIs,2)); % Initialize segmented array
AROI = zeros (1,6); % Array to hold particle areas
binNuc = segs; % Array to hold pre-segmented binarized nuclei for use in vimentin stains

for i = 1:6 % Iterate through each subROI
    %s = graythresh(ROIs(:,:,i));
    %binary_mask = imbinarize(ROIs(:,:,i),'adaptive', 'Sensitivity', s); %
    %Binarize ROI based on calculated threshold
    binary_mask = imbinarize(ROIs(:,:,i)); % Binarize ROI basde on default algorithm
    binary_label = bwlabel(binary_mask, 4); % Label each identified structure in the binary mask
    for labels = 1:max(binary_label,[],'all') % Remove all particles less than .00005 * size 
        count = numel(find(binary_label == labels));
        if (count < 0.00005*numel(ROIs(:,:,i)))
            binary_label(binary_label == labels) = 0;
        end
    end
    binary_mask = (binary_label ~=0);
    
    % Check to set additional particle removal for noisy images - will continue looping until a value is specified 
    binary_areaopen = binary_mask;
    if i == 1 && strcmp(regName, 'IF') % Only perform for 1st image of IF region
        figure; imshowpair(subR1, binary_mask, 'montage'); % Compare grayscale and binary images to determine if area opening is needed
        prompt = {'Remove small particles? Click cancel to skip.'};
        val = str2double(inputdlg(prompt,'Input pixel threshold for removal',[1 35]));
        areaVal = 0; % Initialize variable for px threshold to be used in RM once defined from IF run
        
        % Continue to re-prompt until user exits out of loop manually using cancel - 
        while val ~= 0 % Exit out of loop only by selecting cancel, not entering 0
            areaVal = val; % Saves value for use in RM since val will be overwritten
            binary_areaopen = bwareaopen(binary_mask, val);
            close all;
            montage({subR1, binary_mask, binary_areaopen}, 'Size', [1 3]); % New comparison and prompt
            prompt = {'Try a different threshold? Click cancel to use specified value'};
            val = str2double(inputdlg(prompt,'Input pixel threshold for removal',[1 35]));
        end
    else
        binary_areaopen = bwareaopen(binary_mask, areaVal); % Automatically applies selected threshold to all other subROIs in the region
    end
    binary_mask = binary_areaopen; 
    
    % Calculate area of binarized image prior to segmentation - watershed
    % segmentation can reduce particle area depending on input quality
    AROI(1,i) = sum(binary_mask, 'all'); 
    
    % Perform watershed segmentation
    D = -bwdist(~binary_mask);
    D(~binary_mask) = -Inf;
    L = watershed(D);
    segmentation = imclearborder(L);
    %imshow(label2rgb(segmentation,'jet','k'));
    watershed_label = bwlabel(segmentation,4);
    imshow(watershed_label);
    
    % Add pre- and post-segmentation images to their own arrays
    if i == 1
        segs = cat(3, watershed_label);
        binNuc = cat(3, binary_mask);
    else
        segs = cat(3,segs, watershed_label);
        binNuc = cat(3, binNuc, binary_mask);
    end
end
close all;

%% Calculation of nuclear count and area

% regionprops stats for each subROI
stats1 = regionprops(segs(:,:,1),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');
stats2 = regionprops(segs(:,:,2),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');
stats3 = regionprops(segs(:,:,3),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');
stats4 = regionprops(segs(:,:,4),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');
stats5 = regionprops(segs(:,:,5),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');
stats6 = regionprops(segs(:,:,6),'Area','Centroid','MajorAxisLength','MinorAxisLength','Orientation','PixelList', 'Image');

% Calculate count and area summary statistics, with first line counting nuclei and second summing area
C1 = size(stats1,1);
C2 = size(stats2,1);
C3 = size(stats3,1);
C4 = size(stats4,1);
C5 = size(stats5,1);
C6 = size(stats6,1);
cROI = [C1, C2, C3, C4, C5, C6];

% Scale normalization assumes conversion given is px/um, normalized to #/mm^2
scaleCount_Inner = mean(cROI(1:3)) / (subROI_px/((scale.^2)*1000.^2));
scaleCount_Outer = mean(cROI(4:6)) / (subROI_px/((scale.^2)*1000.^2));
scaleCount = mean(cROI) / (subROI_px/((scale.^2)*1000.^2)); 

A1 = AROI(1) / subROI_px;
A2 = AROI(2) / subROI_px;
A3 = AROI(3) / subROI_px;
A4 = AROI(4) / subROI_px;
A5 = AROI(5) / subROI_px;
A6 = AROI(6) / subROI_px;

% Normalized by analyzed area: Area (px^2)/(px^2)
scaleArea_Inner = mean(AROI(1:3)) / subROI_px;
scaleArea_Outer = mean(AROI(4:6)) / subROI_px;
scaleArea = mean(AROI(1,1:6)) / subROI_px; 

%% Nearest Neighbors Analysis - see nearestNeighbors.m

stats1 = nearestNeighbors(stats1, path, segs(:,:,1), 20, scale, strcat(regName,'_NN1'));
stats2 = nearestNeighbors(stats2, path, segs(:,:,2), 20, scale, strcat(regName,'_NN2'));
stats3 = nearestNeighbors(stats3, path, segs(:,:,3), 20, scale, strcat(regName,'_NN3'));
stats4 = nearestNeighbors(stats4, path, segs(:,:,4), 20, scale, strcat(regName,'_NN4'));
stats5 = nearestNeighbors(stats5, path, segs(:,:,5), 20, scale, strcat(regName,'_NN5'));
stats6 = nearestNeighbors(stats6, path, segs(:,:,6), 20, scale, strcat(regName,'_NN6'));

% Calculations for summary data
% Average number of nuclei within specified search distance
nnC1 = mean([stats1.nnCount]);
nnC2 = mean([stats2.nnCount]);
nnC3 = mean([stats3.nnCount]);
nnC4 = mean([stats4.nnCount]);
nnC5 = mean([stats5.nnCount]);
nnC6 = mean([stats6.nnCount]);
nnCount = [nnC1 nnC2 nnC3 nnC4 nnC5 nnC6];

% Average closest distance of nuclei
nnD1 = mean([stats1.nearestDist_um], 'omitnan');
nnD2 = mean([stats2.nearestDist_um], 'omitnan');
nnD3 = mean([stats3.nearestDist_um], 'omitnan');
nnD4 = mean([stats4.nearestDist_um], 'omitnan');
nnD5 = mean([stats5.nearestDist_um], 'omitnan');
nnD6 = mean([stats6.nearestDist_um], 'omitnan');
nearestDist = [nnD1 nnD2 nnD3 nnD4 nnD5 nnD6];

avg_nnCount = mean(nnCount);
nnCount_Inner = mean(nnCount(1:3));
nnCount_Outer = mean(nnCount(4:6));
avg_nearestDist = mean(nearestDist);
nearestDist_Inner = mean(nearestDist(1:3));
nearestDist_Outer = mean(nearestDist(4:6));

percentNeighbors = [sum(~isnan([stats1.nearestDist_um]))/length([stats1.nearestDist_um]) sum(~isnan([stats2.nearestDist_um]))/length([stats2.nearestDist_um]) sum(~isnan([stats3.nearestDist_um]))/length([stats3.nearestDist_um]) sum(~isnan([stats4.nearestDist_um]))/length([stats4.nearestDist_um]) sum(~isnan([stats5.nearestDist_um]))/length([stats5.nearestDist_um]) sum(~isnan([stats6.nearestDist_um]))/length([stats6.nearestDist_um])];
percentLone = 1 - mean(percentNeighbors); % Number of cells with no nearest neighbors

%% Process Histograms - see histProcess.m

[stats1] = histProcess (region, stats1, ROIpos(1, :), radRing, 1.1, 5, segs(:,:,1), regName, '1', path);
[stats2] = histProcess (region, stats2, ROIpos(2, :), radRing, 1.1, 5, segs(:,:,2), regName, '2', path);
[stats3] = histProcess (region, stats3, ROIpos(3, :), radRing, 1.1, 5, segs(:,:,3), regName, '3', path);
[stats4] = histProcess (region, stats4, ROIpos(4, :), radRing, 1.1, 5, segs(:,:,4), regName, '4', path);
[stats5] = histProcess (region, stats5, ROIpos(5, :), radRing, 1.1, 5, segs(:,:,5), regName, '5', path);
[stats6] = histProcess (region, stats6, ROIpos(6, :), radRing, 1.1, 5, segs(:,:,6), regName, '6', path);

% Calculations for summary data
avgorHist = cat(1, stats1.orHist, stats2.orHist, stats3.orHist, stats4.orHist, stats5.orHist, stats6.orHist);
orHist_Inner = cat(1, stats1.orHist, stats2.orHist, stats3.orHist);
orHist_Outer = cat(1, stats4.orHist, stats5.orHist, stats6.orHist);

% Omit NaN values prior to histogram creation
avgorHist = avgorHist(~isnan(avgorHist)); 
orHist_Inner = orHist_Inner(~isnan(orHist_Inner)); 
orHist_Outer = orHist_Outer(~isnan(orHist_Outer)); 
 
% Save histograms
iptsetpref('ImshowBorder','tight');
histogram(avgorHist, 18, 'Normalization', 'probability');
saveas(gcf,fullfile(path, strcat(regName,'_avgorHist.tif')),'tiff');
histogram(orHist_Inner, 18, 'Normalization', 'probability');
saveas(gcf,fullfile(path, strcat(regName,'_avgorHist_Inner.tif')),'tiff');
histogram(orHist_Outer, 18, 'Normalization', 'probability');
saveas(gcf,fullfile(path, strcat(regName,'_avgorHist_Outer.tif')),'tiff');
close all;

% Output values to QuantHist spreadsheet for compiled histogram creation
if strcmp(regName, 'IF')
    writematrix(avgorHist,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'A:A');
    writematrix(orHist_Inner,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'B:B');
    writematrix(orHist_Outer,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'C:C');
else
    writematrix(avgorHist,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'D:D');
    writematrix(orHist_Inner,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'E:E');
    writematrix(orHist_Outer,'ARA_QuantHist.xlsx', 'Sheet', filename(1:12), 'Range', 'F:F');
end

%% Output data

summary = [scaleCount, scaleCount_Inner, scaleCount_Outer, C1, C2, C3, C4, C5, C6, 0;
    scaleArea, scaleArea_Inner, scaleArea_Outer, A1, A2, A3, A4, A5, A6, 0;
    avg_nnCount, nnCount_Inner, nnCount_Outer, percentLone, nnC1, nnC2, nnC3, nnC4, nnC5, nnC6;
    avg_nearestDist, nearestDist_Inner, nearestDist_Outer, nnD1 nnD2 nnD3 nnD4 nnD5 nnD6 0];

end % End of analyzeNuclei.m
